eReefs Xarray interactiveΒΆ

import os
import numpy as np
import xarray as xr

import cartopy
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
cartopy.config['data_dir'] = os.getenv('CARTOPY_DIR', cartopy.config.get('data_dir'))

import cmocean

import holoviews as hv
from holoviews import opts, dim

import geoviews as gv
import geoviews.feature as gf
from cartopy import crs

gv.extension('bokeh')

import warnings
warnings.filterwarnings("ignore", category=RuntimeWarning) 
year = 2018

# GBR4
base_url = "http://thredds.ereefs.aims.gov.au/thredds/dodsC/s3://aims-ereefs-public-prod/derived/ncaggregate/ereefs/GBR4_H2p0_B3p1_Cq3b_Dhnd/daily-monthly/EREEFS_AIMS-CSIRO_GBR4_H2p0_B3p1_Cq3b_Dhnd_bgc_daily-monthly-"

biofiles = [f"{base_url}{year}-{month:02}.nc" for month in range(1, 2)]
biofiles
['http://thredds.ereefs.aims.gov.au/thredds/dodsC/s3://aims-ereefs-public-prod/derived/ncaggregate/ereefs/GBR4_H2p0_B3p1_Cq3b_Dhnd/daily-monthly/EREEFS_AIMS-CSIRO_GBR4_H2p0_B3p1_Cq3b_Dhnd_bgc_daily-monthly-2018-01.nc']
ds_bio = xr.open_mfdataset(biofiles)
ds_bio
<xarray.Dataset>
Dimensions:          (k: 17, latitude: 723, longitude: 491, time: 31)
Coordinates:
    zc               (k) float64 dask.array<chunksize=(17,), meta=np.ndarray>
  * time             (time) datetime64[ns] 2018-01-01T02:00:00 ... 2018-01-31...
  * latitude         (latitude) float64 -28.7 -28.67 -28.64 ... -7.066 -7.036
  * longitude        (longitude) float64 142.2 142.2 142.2 ... 156.8 156.8 156.9
Dimensions without coordinates: k
Data variables: (12/101)
    alk              (time, k, latitude, longitude) float32 dask.array<chunksize=(31, 17, 723, 491), meta=np.ndarray>
    BOD              (time, k, latitude, longitude) float32 dask.array<chunksize=(31, 17, 723, 491), meta=np.ndarray>
    Chl_a_sum        (time, k, latitude, longitude) float32 dask.array<chunksize=(31, 17, 723, 491), meta=np.ndarray>
    CO32             (time, k, latitude, longitude) float32 dask.array<chunksize=(31, 17, 723, 491), meta=np.ndarray>
    DIC              (time, k, latitude, longitude) float32 dask.array<chunksize=(31, 17, 723, 491), meta=np.ndarray>
    DIN              (time, k, latitude, longitude) float32 dask.array<chunksize=(31, 17, 723, 491), meta=np.ndarray>
    ...               ...
    SGH_N            (time, latitude, longitude) float32 dask.array<chunksize=(31, 723, 491), meta=np.ndarray>
    SGH_N_pr         (time, latitude, longitude) float32 dask.array<chunksize=(31, 723, 491), meta=np.ndarray>
    SGHROOT_N        (time, latitude, longitude) float32 dask.array<chunksize=(31, 723, 491), meta=np.ndarray>
    SGROOT_N         (time, latitude, longitude) float32 dask.array<chunksize=(31, 723, 491), meta=np.ndarray>
    TSSM             (time, latitude, longitude) float32 dask.array<chunksize=(31, 723, 491), meta=np.ndarray>
    Zenith2D         (time, latitude, longitude) float32 dask.array<chunksize=(31, 723, 491), meta=np.ndarray>
Attributes: (12/20)
    Conventions:                     CF-1.0
    NCO:                             netCDF Operators version 4.7.7 (Homepage...
    RunID:                           2
    _CoordSysBuilder:                ucar.nc2.dataset.conv.CF1Convention
    aims_ncaggregate_buildDate:      2020-08-21T23:07:30+10:00
    aims_ncaggregate_datasetId:      products__ncaggregate__ereefs__GBR4_H2p0...
    ...                              ...
    paramfile:                       /home/bai155/EMS_solar2/gbr4_H2p0_B3p1_C...
    paramhead:                       eReefs 4 km grid. SOURCE Catchments with...
    technical_guide_link:            https://eatlas.org.au/pydio/public/aims-...
    technical_guide_publish_date:    2020-08-18
    title:                           eReefs AIMS-CSIRO GBR4 BioGeoChemical 3....
    DODS_EXTRA.Unlimited_Dimension:  time
base_url2 = "http://thredds.ereefs.aims.gov.au/thredds/dodsC/s3://aims-ereefs-public-prod/derived/ncaggregate/ereefs/gbr4_v2/daily-monthly/EREEFS_AIMS-CSIRO_gbr4_v2_hydro_daily-monthly-"
hydrofiles = [f"{base_url2}{year}-{month:02}.nc" for month in range(1, 2)]
hydrofiles
['http://thredds.ereefs.aims.gov.au/thredds/dodsC/s3://aims-ereefs-public-prod/derived/ncaggregate/ereefs/gbr4_v2/daily-monthly/EREEFS_AIMS-CSIRO_gbr4_v2_hydro_daily-monthly-2018-01.nc']
ds_hydro = xr.open_mfdataset(hydrofiles)
ds_hydro
<xarray.Dataset>
Dimensions:      (k: 17, latitude: 723, longitude: 491, time: 31)
Coordinates:
    zc           (k) float64 dask.array<chunksize=(17,), meta=np.ndarray>
  * time         (time) datetime64[ns] 2017-12-31T14:00:00 ... 2018-01-30T14:...
  * latitude     (latitude) float64 -28.7 -28.67 -28.64 ... -7.096 -7.066 -7.036
  * longitude    (longitude) float64 142.2 142.2 142.2 ... 156.8 156.8 156.9
Dimensions without coordinates: k
Data variables:
    mean_cur     (time, k, latitude, longitude) float32 dask.array<chunksize=(31, 17, 723, 491), meta=np.ndarray>
    salt         (time, k, latitude, longitude) float32 dask.array<chunksize=(31, 17, 723, 491), meta=np.ndarray>
    temp         (time, k, latitude, longitude) float32 dask.array<chunksize=(31, 17, 723, 491), meta=np.ndarray>
    u            (time, k, latitude, longitude) float32 dask.array<chunksize=(31, 17, 723, 491), meta=np.ndarray>
    v            (time, k, latitude, longitude) float32 dask.array<chunksize=(31, 17, 723, 491), meta=np.ndarray>
    mean_wspeed  (time, latitude, longitude) float32 dask.array<chunksize=(31, 723, 491), meta=np.ndarray>
    eta          (time, latitude, longitude) float32 dask.array<chunksize=(31, 723, 491), meta=np.ndarray>
    wspeed_u     (time, latitude, longitude) float32 dask.array<chunksize=(31, 723, 491), meta=np.ndarray>
    wspeed_v     (time, latitude, longitude) float32 dask.array<chunksize=(31, 723, 491), meta=np.ndarray>
Attributes: (12/19)
    Conventions:                     CF-1.0
    Run_ID:                          2
    _CoordSysBuilder:                ucar.nc2.dataset.conv.CF1Convention
    aims_ncaggregate_buildDate:      2020-08-21T12:50:07+10:00
    aims_ncaggregate_datasetId:      products__ncaggregate__ereefs__gbr4_v2__...
    aims_ncaggregate_firstDate:      2018-01-01T00:00:00+10:00
    ...                              ...
    paramhead:                       GBR 4km resolution grid
    shoc_version:                    v1.1 rev(5620)
    technical_guide_link:            https://eatlas.org.au/pydio/public/aims-...
    technical_guide_publish_date:    2020-08-18
    title:                           eReefs AIMS-CSIRO GBR4 Hydrodynamic v2 d...
    DODS_EXTRA.Unlimited_Dimension:  time
reef_lat = -18.82
reef_lon = 147.64
min_lon = 146
min_lat = -21
max_lon = 150
max_lat = -16

lon_bnds = [min_lon, max_lon]
lat_bnds = [min_lat, max_lat]
ds_bio_clip = ds_bio.sel(latitude=slice(*lat_bnds), longitude=slice(*lon_bnds))
ds_hydro_clip = ds_hydro.sel(latitude=slice(*lat_bnds), longitude=slice(*lon_bnds))
zcvar = -1
new_ds = ds_hydro_clip.isel(k=zcvar)
new_ds
<xarray.Dataset>
Dimensions:      (latitude: 167, longitude: 134, time: 31)
Coordinates:
    zc           float64 dask.array<chunksize=(), meta=np.ndarray>
  * time         (time) datetime64[ns] 2017-12-31T14:00:00 ... 2018-01-30T14:...
  * latitude     (latitude) float64 -20.99 -20.96 -20.93 ... -16.04 -16.01
  * longitude    (longitude) float64 146.0 146.0 146.1 ... 149.9 150.0 150.0
Data variables:
    mean_cur     (time, latitude, longitude) float32 dask.array<chunksize=(31, 167, 134), meta=np.ndarray>
    salt         (time, latitude, longitude) float32 dask.array<chunksize=(31, 167, 134), meta=np.ndarray>
    temp         (time, latitude, longitude) float32 dask.array<chunksize=(31, 167, 134), meta=np.ndarray>
    u            (time, latitude, longitude) float32 dask.array<chunksize=(31, 167, 134), meta=np.ndarray>
    v            (time, latitude, longitude) float32 dask.array<chunksize=(31, 167, 134), meta=np.ndarray>
    mean_wspeed  (time, latitude, longitude) float32 dask.array<chunksize=(31, 167, 134), meta=np.ndarray>
    eta          (time, latitude, longitude) float32 dask.array<chunksize=(31, 167, 134), meta=np.ndarray>
    wspeed_u     (time, latitude, longitude) float32 dask.array<chunksize=(31, 167, 134), meta=np.ndarray>
    wspeed_v     (time, latitude, longitude) float32 dask.array<chunksize=(31, 167, 134), meta=np.ndarray>
Attributes: (12/19)
    Conventions:                     CF-1.0
    Run_ID:                          2
    _CoordSysBuilder:                ucar.nc2.dataset.conv.CF1Convention
    aims_ncaggregate_buildDate:      2020-08-21T12:50:07+10:00
    aims_ncaggregate_datasetId:      products__ncaggregate__ereefs__gbr4_v2__...
    aims_ncaggregate_firstDate:      2018-01-01T00:00:00+10:00
    ...                              ...
    paramhead:                       GBR 4km resolution grid
    shoc_version:                    v1.1 rev(5620)
    technical_guide_link:            https://eatlas.org.au/pydio/public/aims-...
    technical_guide_publish_date:    2020-08-18
    title:                           eReefs AIMS-CSIRO GBR4 Hydrodynamic v2 d...
    DODS_EXTRA.Unlimited_Dimension:  time
week_ds = new_ds.resample(time='1W').mean(dim='time').drop('zc')[['mean_cur','u','v']]
week_ds
<xarray.Dataset>
Dimensions:    (latitude: 167, longitude: 134, time: 6)
Coordinates:
  * time       (time) datetime64[ns] 2017-12-31 2018-01-07 ... 2018-02-04
  * latitude   (latitude) float64 -20.99 -20.96 -20.93 ... -16.07 -16.04 -16.01
  * longitude  (longitude) float64 146.0 146.0 146.1 146.1 ... 149.9 150.0 150.0
Data variables:
    mean_cur   (time, latitude, longitude) float32 dask.array<chunksize=(1, 167, 134), meta=np.ndarray>
    u          (time, latitude, longitude) float32 dask.array<chunksize=(1, 167, 134), meta=np.ndarray>
    v          (time, latitude, longitude) float32 dask.array<chunksize=(1, 167, 134), meta=np.ndarray>
load_week = week_ds.load()
load_week
---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-11-9de0457538f4> in <module>
      1 load_week = week_ds.load()
----> 2 load_week

/usr/share/miniconda/envs/envireef/lib/python3.8/site-packages/IPython/core/displayhook.py in __call__(self, result)
    260             self.start_displayhook()
    261             self.write_output_prompt()
--> 262             format_dict, md_dict = self.compute_format_data(result)
    263             self.update_user_ns(result)
    264             self.fill_exec_result(result)

/usr/share/miniconda/envs/envireef/lib/python3.8/site-packages/IPython/core/displayhook.py in compute_format_data(self, result)
    149 
    150         """
--> 151         return self.shell.display_formatter.format(result)
    152 
    153     # This can be set to True by the write_output_prompt method in a subclass

/usr/share/miniconda/envs/envireef/lib/python3.8/site-packages/IPython/core/formatters.py in format(self, obj, include, exclude)
    178             md = None
    179             try:
--> 180                 data = formatter(obj)
    181             except:
    182                 # FIXME: log the exception

<decorator-gen-2> in __call__(self, obj)

/usr/share/miniconda/envs/envireef/lib/python3.8/site-packages/IPython/core/formatters.py in catch_format_error(method, self, *args, **kwargs)
    222     """show traceback on failed format call"""
    223     try:
--> 224         r = method(self, *args, **kwargs)
    225     except NotImplementedError:
    226         # don't warn on NotImplementedErrors

/usr/share/miniconda/envs/envireef/lib/python3.8/site-packages/IPython/core/formatters.py in __call__(self, obj)
    343             method = get_real_method(obj, self.print_method)
    344             if method is not None:
--> 345                 return method()
    346             return None
    347         else:

/usr/share/miniconda/envs/envireef/lib/python3.8/site-packages/xarray/core/dataset.py in _repr_html_(self)
   1808         if OPTIONS["display_style"] == "text":
   1809             return f"<pre>{escape(repr(self))}</pre>"
-> 1810         return formatting_html.dataset_repr(self)
   1811 
   1812     def info(self, buf=None) -> None:

/usr/share/miniconda/envs/envireef/lib/python3.8/site-packages/xarray/core/formatting_html.py in dataset_repr(ds)
    281     sections = [
    282         dim_section(ds),
--> 283         coord_section(ds.coords),
    284         datavar_section(ds.data_vars),
    285         attr_section(ds.attrs),

/usr/share/miniconda/envs/envireef/lib/python3.8/site-packages/xarray/core/formatting_html.py in _mapping_section(mapping, name, details_func, max_items_collapse, enabled)
    171     return collapsible_section(
    172         name,
--> 173         details=details_func(mapping),
    174         n_items=n_items,
    175         enabled=enabled,

/usr/share/miniconda/envs/envireef/lib/python3.8/site-packages/xarray/core/formatting_html.py in summarize_coords(variables)
     90     coords = {}
     91     for k, v in variables.items():
---> 92         coords.update(**summarize_coord(k, v))
     93 
     94     vars_li = "".join(f"<li class='xr-var-item'>{v}</li>" for v in coords.values())

/usr/share/miniconda/envs/envireef/lib/python3.8/site-packages/xarray/core/formatting_html.py in summarize_coord(name, var)
     84             return coords
     85 
---> 86     return {name: summarize_variable(name, var, is_index)}
     87 
     88 

/usr/share/miniconda/envs/envireef/lib/python3.8/site-packages/xarray/core/formatting_html.py in summarize_variable(name, var, is_index, dtype, preview)
    112     preview = preview or escape(inline_variable_array_repr(variable, 35))
    113     attrs_ul = summarize_attrs(var.attrs)
--> 114     data_repr = short_data_repr_html(variable)
    115 
    116     attrs_icon = _icon("icon-file-text2")

/usr/share/miniconda/envs/envireef/lib/python3.8/site-packages/xarray/core/formatting_html.py in short_data_repr_html(array)
     26         return internal_data._repr_html_()
     27     else:
---> 28         text = escape(short_data_repr(array))
     29         return f"<pre>{text}</pre>"
     30 

/usr/share/miniconda/envs/envireef/lib/python3.8/site-packages/xarray/core/formatting.py in short_data_repr(array)
    481         return limit_lines(repr(array.data), limit=40)
    482     elif array._in_memory or array.size < 1e5:
--> 483         return short_numpy_repr(array)
    484     else:
    485         # internal xarray array type

/usr/share/miniconda/envs/envireef/lib/python3.8/site-packages/xarray/core/formatting.py in short_numpy_repr(array)
    470     options["edgeitems"] = edgeitems
    471     with set_numpy_options(**options):
--> 472         return repr(array)
    473 
    474 

/usr/share/miniconda/envs/envireef/lib/python3.8/site-packages/numpy/core/arrayprint.py in _array_repr_implementation(arr, max_line_width, precision, suppress_small, array2string)
   1411         lst = repr(arr.item())
   1412     elif arr.size > 0 or arr.shape == (0,):
-> 1413         lst = array2string(arr, max_line_width, precision, suppress_small,
   1414                            ', ', prefix, suffix=suffix)
   1415     else:  # show zero-length shape unless it is (0,)

/usr/share/miniconda/envs/envireef/lib/python3.8/site-packages/numpy/core/arrayprint.py in array2string(a, max_line_width, precision, suppress_small, separator, prefix, style, formatter, threshold, edgeitems, sign, floatmode, suffix, legacy)
    690         return "[]"
    691 
--> 692     return _array2string(a, options, separator, prefix)
    693 
    694 

/usr/share/miniconda/envs/envireef/lib/python3.8/site-packages/numpy/core/arrayprint.py in wrapper(self, *args, **kwargs)
    466             repr_running.add(key)
    467             try:
--> 468                 return f(self, *args, **kwargs)
    469             finally:
    470                 repr_running.discard(key)

/usr/share/miniconda/envs/envireef/lib/python3.8/site-packages/numpy/core/arrayprint.py in _array2string(a, options, separator, prefix)
    499     next_line_prefix += " "*len(prefix)
    500 
--> 501     lst = _formatArray(a, format_function, options['linewidth'],
    502                        next_line_prefix, separator, options['edgeitems'],
    503                        summary_insert, options['legacy'])

/usr/share/miniconda/envs/envireef/lib/python3.8/site-packages/numpy/core/arrayprint.py in _formatArray(a, format_function, line_width, next_line_prefix, separator, edge_items, summary_insert, legacy)
    843     try:
    844         # invoke the recursive part with an initial index and prefix
--> 845         return recurser(index=(),
    846                         hanging_indent=next_line_prefix,
    847                         curr_width=line_width)

/usr/share/miniconda/envs/envireef/lib/python3.8/site-packages/numpy/core/arrayprint.py in recurser(index, hanging_indent, curr_width)
    799 
    800             for i in range(trailing_items, 1, -1):
--> 801                 word = recurser(index + (-i,), next_hanging_indent, next_width)
    802                 s, line = _extendLine_pretty(
    803                     s, line, word, elem_width, hanging_indent, legacy)

/usr/share/miniconda/envs/envireef/lib/python3.8/site-packages/numpy/core/arrayprint.py in recurser(index, hanging_indent, curr_width)
    753 
    754         if axes_left == 0:
--> 755             return format_function(a[index])
    756 
    757         # when recursing, add a space to align with the [ added, and reduce the

/usr/share/miniconda/envs/envireef/lib/python3.8/site-packages/numpy/core/arrayprint.py in __call__(self, x)
    992                                       exp_digits=self.exp_size)
    993         else:
--> 994             return dragon4_positional(x,
    995                                       precision=self.precision,
    996                                       unique=self.unique,

KeyboardInterrupt: 
lat = load_week.latitude
lon = load_week.longitude

# Convert to magnitude and angle
mag = np.sqrt(load_week.u**2 + load_week.v**2)
angle = (np.pi/2.) - np.arctan2(load_week.u/mag, load_week.v/mag)
load_week["mag"] = (['time', 'latitude', 'longitude'],  mag)
load_week["angle"] = (['time', 'latitude', 'longitude'],  angle)
dataset = gv.Dataset(load_week, ['longitude', 'latitude', 'time'], 
                     'mag', crs=crs.PlateCarree())
images = dataset.to(gv.Image,dynamic=True)
coastline = gf.coastline(line_width=2,line_color='k').opts(projection=crs.PlateCarree(),scale='10m')
land = gf.land.options(scale='10m', fill_color='lightgray')
resample = load_week.isel(longitude=slice(None, None, 7), latitude=slice(None, None, 7))
resample
<xarray.Dataset>
Dimensions:    (latitude: 24, longitude: 20, time: 6)
Coordinates:
  * time       (time) datetime64[ns] 2017-12-31 2018-01-07 ... 2018-02-04
  * latitude   (latitude) float64 -20.99 -20.78 -20.57 ... -16.58 -16.37 -16.16
  * longitude  (longitude) float64 146.0 146.2 146.4 146.6 ... 149.6 149.8 150.0
Data variables:
    mean_cur   (time, latitude, longitude) float32 nan nan nan ... 0.2521 0.1418
    u          (time, latitude, longitude) float32 nan nan ... -0.0868 -0.0706
    v          (time, latitude, longitude) float32 nan nan ... -0.07677 0.004869
    mag        (time, latitude, longitude) float32 nan nan ... 0.1159 0.07077
    angle      (time, latitude, longitude) float32 nan nan nan ... 3.866 3.073
resample.mag.max()
<xarray.DataArray 'mag' ()>
array(0.98175, dtype=float32)
lat = resample.latitude
lon = resample.longitude
# label = 'Arrows scale with plot width, not view'
# vectorfield = hv.VectorField((lon, lat, resample.angle[0,:,:], resample.mag[0,:,:]))
# vectorfield.relabel(label)
# kdims=[hv.Dimension('step', values=np.arange(0,6)),
#        hv.Dimension('time', values=resample.time.values)]
# hmap = hv.HoloMap(vfield, kdims=kdims)
#hmap = hv.DynamicMap(vfield, kdims=kdims)
max_mag = resample.mag.max()
nb = resample.angle.shape[0]

a_var = resample.time.values
y_list = []
for i in range(nb):
    y_list.append(gv.VectorField((lon, lat, resample.angle[i,:,:], resample.mag[i,:,:]/max_mag), crs=crs.PlateCarree()))

# create HoloMap object    
dict_y = {a_var[i]:y_list[i] for i in range(nb)}

hmap = hv.HoloMap(dict_y, kdims="time").opts(opts.VectorField(magnitude=dim('Magnitude')*0.25, color='k', width=600, height=500,
                                                              pivot='tip', line_width=1, title='eReefs GBR 4km', 
                                                              rescale_lengths=False,
                                                              projection=crs.PlateCarree()))

# hmap
images.opts(cmap=cmocean.cm.speed, colorbar=True, width=500, height=500, clim=(0.1,1.0)) * coastline * land * hmap
figs = images.opts(cmap=cmocean.cm.speed, colorbar=True, width=500, height=500) * coastline * land * hmap
hv.save(figs, 'graph.html')
                                             
# GOTO: https://raw.githack.com
# https://github.com/TristanSalles/EnviReef/main/book/xarray/graph.html
from IPython.display import IFrame
from IPython.core.display import display
display(IFrame('https://raw.githack.com/TristanSalles/EnviReef/main/book/xarray/graph.html', width='950px', height='500px'))